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I.  INTRODUCTION 


Among  the  various  methods  for  the  representation  of  the 
earth's  anomalous  gravity  field,  the  point  mass  model 
representation  is  particularly  attractive  both  because  of 
its  conceptual  simplicity  and  probably  because  of  its  closest 
"neighborhood"  to  geological  reality  among  all  proposed  methods. 
The  generation  of  a  point  mass  model,  however,  is  not  only 
non-unique,  but  also  not  quite  simple.  The  difficulty  lies 
in  a)  the  non-uniqueness  cf  the  problem,  b)  the  computa¬ 
tionally  demanding  relation  between  mean  values  and  point 

values,  c)  the  use  of  mean  values  of  various  kinds  like  5°  x  5°, 

.0  ,o  . 

lxl,  and  so  on. 

In  a  former  report  (Siinkel,  1982)  we  demonstrated  the  reasons 
for  a  multi-level  mass  point  model  employing  known  statistical 
properties  of  the  earth's  gravity  field.  In  the  present  contri¬ 
bution  the  conceptual  mathematical  problems  have  been  solved 
and  a  conceivable  algorithm  for  the  actual  generation  of  a 
mass  point  model  presented.  Frequency  domain  methods  have 
been  frequently  used  throughout  this  paper. 

The  evaluation  of  the  kernel  which  relates  point  masses  and 
mean  gravity  disturbances,  requires  integration  on  the  sphere 
over  a  limited  area  -  a  time  -  consuming  process.  In  order 
to  speed  up  calculation,  a  fast  approximation  has  to  be 
designed  and  the  approximation  error  to  be  estimated.  Using 
Peano's  theorem  on  the  sphere  and  the  method  of  Sard,  "best" 
approximations  have  been  derived  for  various  kernel  approxi¬ 
mation  functions  and  the  corresponding  approximation  errors 
have  been  estimated. 

The  relation  between  the  depth  of  the  point  mass  level  and 
a  set  of  mean  gravity  data  at  zero  level  has  been  derived 
from  the  principle  of  minimum  deviation  of  2  operators: 
the  smoothing  operator  (which  transforms  point  data  into 
mean  values)  and  the  operator  which  turns  point  masses  at 


a  certain  level  into  its  gravitation. 


The  algorithm  for  the  mass  point  model  generation  is  entirely 
based  on  fast  Fourier  transform  methods: 

the  data  and  the  operator  are  transformed  into  the  frequency 
domain,  the  solution  of  the  linear  system  performed,  and 
the  result  re- trans formed  into  the  space  domain  yielding 
the  required  set  of  point  masses.  Due  to  the  use  of  mean 
values  of  various  kind,  a  recursive  procedure  had  to  be 
desi gned. 
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2.  POINT  MASSES  AND  MEAN  GRAVITY  DISTURBANCES 


* 


The  gravitational  potential  T  generated  by  a  set  of  point 
masses  {m..}  is  given  by 

j  Gm . 

T(p)  *  ;  (2-1) 

G  denotes  the  gravitational  constant,  1(P,Q  )  the  spatial  dis¬ 
tance  between  the  calculation  point  P  and  the  location  Q ^ 
of  the  point  mass  m^  .  The  negative  radial  component  of  the 
corresponding  gravitational  vector  VT  will  be  denoted  gravity 
disturbance  and  can  easily  be  derived  from  (2-1), 


,T  J  r  -r  . cosvp  . 

«9(P)  ■  -  fl  -  I  -V5 - 1  Gm- 

3r  1=1  13(P,Qj> 


(2-2) 


where  r  stands  for  radius.  Given  {m^}  ,  the  calculation  of 
6g  is  straightforward. 

Let  us  now  investigate  the  relation  between  point  masses 
and  mean  gravity  disturbances,  assuming  that  all  the  point  masses 
are  located  on  a  sphere  with  radius  a  <  1  with  the  gravity  dis¬ 
turbance  determined  for  a  point  on  the  unit  sphere  a  =  1  con¬ 
centric  to  the  former  sphere.  Then  rp  =  1  and  r^  =  a  and  the 
position  of  P  on  the  unit  sphere  is  given  by  the  unit  vector 
5  ,  the  position  of  by  on  with  the  unit  vector  n  . 

The  kernel  in  equation  (2-2)  is  homogeneous  and  isotropic 
and  can  be  represented  in  terms  of  a  series  of  Legendre  polyno¬ 
mials  Pn(cosij/)  (Heiskanen  &  Moritz,  1967,  p.35). 


l  (n  +  1  )anPn  (  COSiJ; ) 
n =0 


(2-3) 
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Representing  Pn  in  terms  of  orthonormal  surface  spherical  har¬ 
monics  { 4>  }  ,  n  =  0,...;  m  =  -n,...,  n, 

no 

I'nJtKlfOdoU)  -  +)  (2-4) 

a 

with  the  decomposition  formula 


P 

n 


U-n) 


“  TnTT  l  *nra(^Un) 


(2-5) 


the  dependence  of  the  kernel  (2-3)  on  the  position  of  P  and 
Q  on  the  unit  sphere,  represented  by  the  unit  vectors  £  and 
n  ,  is  given  by 


(1-qCOSi t  _)  “  _,i  n 

-73 - 1  Sh?T  •“  I 

1  pQ  n=0  tn  =  -n 


(2-6) 


If  we  want  to  determine  the  kernel  between  a  point  mass  and  a 
mean  gravity  disturbance,  we  have  to  calculate  the  mean  of  (2-6) 
with  respect  to  P  , 


([  l-aCOSiJ>  )>  «®  n 

Mp  (-71 - “-}  '  l  &7T  •"  l 

'•1  •  n  =  0  m=-n 


where  M  stands  for  the  mean  value  operator 


M  {(•)}  =  u  *  ( • )  da  (2-8) 

U) 

taken  over  a  limited  part  uca  of  the  unit  sphere  a  .  The  mean 
gravity  disturbance  6g  , 


+)  Note  the  difference  between  our  definition  and  the  definition 
of  fully  normalized  spherical  harmonics  of  (Heiskanen  & 
Moritz,  p.31  ) 


can  then  be  expressed  by 


It  is  evident  that  the  rate  of  convergence  of  (2-9)  depend'  in 
a  and  therefore  on  the  depth  (1-a)  of  the  point  mass  spf  e . 
Shallow  point  masses  will  cause  a  slow  convergence,  deep  >  (nt 
masses  a  rapid  convergence.  (If  a  =  0  ,  all  the  point  ma 
are  concentrated  at  the  center  of  the  sphere  and  all  terms  .1  >  0 
are  anni  h  i  1  a  ted  by  the  powers  of  a  .  )  Therefore,  for  very  deep 
masses  the  summation  can  be  terminated  at  a  modest  n  =  N  ;  how¬ 
ever,  in  our  study  even  for  the  most  favourable  depth  (D  =  550km 
the  summation  has  to  be  carried  out  at  least  up  to  the  degree 
n  a  150  in  order  to  obtain  6  significant  digits.  Consequent¬ 
ly,  the  use  of  equation  (2-9)  for  the  calculation  of  mean  gravity 
disturbances  is  prohibitive  despite  the  existence  of  a  very  smart 
algorithm  for  the  calculation  of  integrals  of  associated  Legen¬ 
dre  functions  (Gerstl  ,  1980). 

As  an  alternative,  we  suggest  a  local  approximation  of  the 
kernel  by  a  low  degree  polynomial  which,  restricted  to  the  sphere 
is  a  1  i near  combi nati on  of  low  degree  spherical  harmonics.  Such 
an  approximation  allows  us  to  estimate  the  error  of  the  kernel's 
mean  value  for  the  area  in  consideration  by  taking  advantage  of 
Peano's  theorem  and  moreover,  to  estimate  the  best  possible 
approximation  in  the  sense  of  Sard. 
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3.  LOCAL  KERNEL  APPROXIMATIONS 


In  a  paper  by  Siinkel  and  Rummel  (1981)  one  possible  local 
kernel  approximation  has  been  proposed;  an  integration  error 
estimate  has  been  given  for  the  very  special  case  that  the  cal¬ 
culation  point  coincides  with  the  pole.  In  this  particular  case, 
the  one-dimensional  Peano  theorem  can  be  applied.  Considerable 
effort  in  overcoming  this  deficiency  has  been  made  by  W.  Freeden 
culminating  in  a  particul arly  beauti  ful  elaborate  investigation 
"On  Spherical  Spline  Interpolation  and  Approximation"  (1981). 

In  the  sequel,  the  results  of  this  paper  are  used  from  an 
application  point  of  view.  In  particular,  we  are  interested  in 
both  the  estimate  of  the  error  which  we  commit,  if  for  the  inte¬ 
gration  a  local  low  degree  polynomial  approximation  of  the  kernel 
is  employed,  and  in  the  best  approximation  by  such  a  polynomial. 
Since  the  author  of  this  paper  favorizes  the  inductive  approach, 
the  simplest  case  will  be  studied  first.  For  the  sake  of  simpli¬ 
city  we  further  assume  in  the  following  derivations  that  the  cal¬ 
culation  point  is  kept  fixed;  then  the  kernel  to  be  studied  is 
a  function  of  the  integration  point  only. 

3.1.  Approximation  of  degree  zero  (J  =  0) 

The  simplest  approximation  of  a  mean  value  of  a  function  f 
is  obviously  given  by  a  single  function  value  which  is  usually 
taken  at  the  center  of  the  area  in  consideration.  The  mean  value 
of  a  function  (for  a  limited  area  w  C  o  of  the  unit  sphere) 
is  a  linear  functional  L  applied  to  the  function  f  and  can 
be  represented  in  terms  of 


.  •  i  •  »  •  , 


,V 


with  the  characteristic  g 


9(n) 


for  r\  e  in 
else  . 


(3-2) 


This  mean  value  is  approximated  by  a  single  function  value  at 
a  preselected  point  ,  which  is  again  a  linear  functional  L 

Lf  =  atjfU^  ;  (3-3) 


the  remainder  Rf  (another  linear  functional)  is  consequently 
given  by 


Rf  -  Lf  -  Lf  -  (L  -  L ) f 


“  Jg(n)f(n)do(n)  -  a x f ( C A )  •  (3-4) 

a 

For  the  following  we  assume  that  feC2(o)  ,  the  space  of  all 
twice  continuously  differentiable  functions  on  a  . 


Peano's  Theorem 


Let  us  assume  that  R  annihilates  all  zero  degree  poly¬ 
nomials,  Rh  =  0  vh  G  IP  .  Then,  according  to  Peano's  theorem 
Rf  can  be  represented  in  terms  of 

Rf  =  [k  (n)A*f(n)do(n)  (3-5) 

I  °  n 
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with  the  spherical  Peano  kernel 

K0(n)  -  -  ^  R5G<'>  (6, n)  (3-6) 

and  Green's  function  (Freeden,  p.566) 


2n  +  l 


C  U,n)  =•  I  jif=5h  P„(6.n)  . 


(3-7) 


In  this  context  A*  denotes  the  Be  1  trami-operator ,  which  is  simp¬ 


ly  the  restriction  of  the  Lapl ace-operator  A  to  the  unit 
sphere  a  , 


n  (1-t2)  te  ♦  77?  [jt]  '  t:  -  cos9  (3'8) 


applied  at  the  point  z  e  a  having  spherical  polar  coordinates 
e  and  1  . 

Using  Schwarz's  inequality,  an  upper  bound  for  the  approxi¬ 
mation  error  can  be  given  by 


(Rf)2  £  j[Ko(n)]2dc(n)  •  J[An*f (n)]2da(n) 


(3-9) 


The  first  integral  on  the  right  hand  side  can  be  shown  (Freeden, 
p.558)  to  be  equal  to 


A:  =  [  K  (  n )  ]  2  da  ( n  )  =  — MG?  U,c) 

°  (4*)2  C  C  ° 


(3-10) 


with  the  iterated  Green's  function 


Gf  U,c)  =  M  2r?t-- - 2  pn(5*c)  • 

°  n=l[n(n+l)]2  n 


(3-11) 


According  to  the  decomposition  theorem,  the  Legendre  polynomials 

P  (C.c)  can  be  expressed  in  terms  of  the  normalized  surface 

spherical  harmonics  $  ,  m  =  -n,...,  n, 

nm 


g?’u,o=  <4,n 


n=  1  [n(n  +  l)p  m=-n 


l  <>  UU  U)  • 

nm  nm 


(3-11 )  * 


The  application  of  (3-10)  yields 


l  - - - 

n= l [ n ( n+1 ) ] 


n 


l  (♦*  “  2*a, ) 

nm  nm  1  nm  1 

m=-n 


(3-10) 1 


The  second  .integral  on  the  right  hand  side  is  equally  simple  to 
evaluate:  since  our  original  kernel  f  to  be  investigated  is 
homogeneous  and  isotropic  and  since  the  integration  in  (3-9) 
has  to  be  carried  out  over  the  whole  unit  sphere  o  ,  the  posi¬ 
tion  of  the  calculation  point  is  immaterial  therefore,  we  choose 
the  pole  9  =  0  .  In  this  case  the  Bel trami -operator  reduces  to 
the  Legendre-operator 


(3-12) 


and  with  our  expression  (2-3)  for  the  kernel  and  with  the  use 
of  the  Legendre  differential  equation 

Zl  d-t2)  Zl  P„(t)  -  -  "(n+l)Pn(t)  (3-13) 

the  operation  A*f  yields 

oo 

A*f  =  -  l  n ( n+1 )2anP  (t)  .  (3-14) 

—  —  i  ^ 


Taking  into  account  the  orthogonality  of  the  Legendre  polynomials 


nm 


(3-15) 


JPn(t)pm(t)dt 


77m  5 


- 1 

(«n  ...  Kronecker  symbol),  we  obtain 

B:=  f  [a*f  (  n)  l2do(  r,)  =  4,  j  (3-16) 

>  n  =  1 

a 

and  the  estimation  for  the  upper  bound  of  the  approximation  error 


|  Rf  |  s  i/AB 


(3-17) 


Let  us  investigate  this  estimate  in  more  detail:  for  the  calcu¬ 
lation  of  A  we  have  a)  to  calculate  the  mean  value  of  surface 
spherical  harmonics  for  the  area  in  consideration  and  b)  to 
evaluate  the  surface  spherical  harmonics  at  a  (single)  point. 

In  this  respect  there  is  practically  no  difference  between  (3-10)' 

—  2 

and  (2-7);  however,  due  to  the  strong  damping  factor  (n(n+l)] 
the  series  (3-10)'  converges  strongly  and,  as  a  consequence,  only 
relatively  few  terms  are  necessary  compared  to  (2-7)  which,  even 
under  optimal  circumstances  (D  =  550km)  ,  has  to  be  evaluated 
up  to  a  very  high  degree  (about  3  times  higher  if  the  result  is 
to  be  accurate  to  6  significant  digits);  the  situation  is  getting 
even  worse  with  smaller  D  .  The  calculation  of  the  factor  B 
is  very  simple  and  has  to  be  carried  out  only  once  for  each  D  . 

It  is  evident  that  the  convergence  of  (3-16)  is  faster  for  a  lar¬ 
ger  D  .  The  behavior  of  (3-16)  is  quite  similar  to  the  one  shown 
in  (Slinkel  ,  1981,  Fig. 2.1,  pp.  9,10)  with  the  maximum  shifted 
to  the  right  due  to  the  higher  power  in  n  . 

The  next  step  will  be  to  find  the  best  value  for  the  coef¬ 
ficient  at  such  that  the  remainder  (3-17)  is  minimized  under 
the  constraint  Rh  =  0  y  h  €  IP 


Since  B  is  independent  of  a]  ,  it  is  a  mere  scale  fac¬ 
tor  and  can  be  put  equal  to  1  .  Therefore,  the  optimization 
problem  can  be  formulated  as  follows: 


_ 3 

9a 


9 

9  A 


-  {7  "  Jg(n)<^0j0(n)do(n)  -  a|*0>0tei) 

1  0 

{ 


}- 


}  =  °  . 


(3-18) 


_  1 

With  q  =  ( 4ir )  2  and  (3-2)  the  first  expression  simplifies  to 

?  -  ‘d-3,)]  *  0  ' 


With  the  constraint  mentioned  before  (second  line  of  (3-18)),  this 


leads  to  the  following  linear  system  with  the  two  unknowns,  a^  , 
and  the  Lagrange  multiplier  A 


l  — r  1  <L  («,)  -  1 


n=*l  [n(n+l)]  m=-n 
1 


l 


1 


L .  rn#_.iv|2  L  nm  run'  I 
n=l  [n(n+l)l  m=-n 

1 


(3-19) 


Since 


n 


l  ♦„»(')  = 

m«-n 


2n+l 


T7 


ve  e  o 


and 


l 


2n+l 


n= 1  [n(n+l)] 


=  »!,  [7  '  ITT?]  =  1 


(3-20) 


the  above  1  inear  system  reduces  to 


1 


1 


4  rr 


Li  .  oj  LO 


with  the  solution 


n=i  [n(n+l  )]2 


l  Wn£l> 


1 


(3-21) 


CO 


+  l 

n  =  1 


1 

[ n ( n  +  1 )  ]  2 


n 


(3-22) 


Therefore,  the  best  approximation  of  the  kernels  integral  by  a 
single  kernel  function  value  is  obtained  by 


Lf  »  f(5,) 


and  the  a  posteriori  error  estimate  of  this  optimal  solution  is 
given  by 


(Rf)2 


00  .  n 

<;/U  l  - I - 2-  l 

n=  1  [  n(n+l)]  ra=-n 


(3-23) 


A  special  case  is  obviously  given  if  the  area  of  averaging  u> 

reduces  to  zero  and  if  c,  is  (as  usual)  located  at  the  center 

of  the  area.  In  this  case  5  =  <J>  (5.)  and  the  error  becomes 

n  m.  nm  i 

evidently  zero  as  one  should  expect.  This  concludes  the  error 
estimation  of  a  zero  degree  local  kernel  approximation. 


How  do  our  formulas  change  if  Lf  is  approximated  by  a 
linear  combination  of  I  >  1  function  values. 


(3-24) 


If  =  l  a±f (C±) 

i=  1  1  1 

such  that  Rh  =  ( L - L ) h  =0  vh  €  Pq  ?  As  a  matter  of  fact  the 
definition  of  the  remainder  will  now  be 


Rf  =  fg(n)f (n)do(n)  -  f  a.f(e  ) 

J  i=  i 


(3-25) 


and  as  a  consequence 

oo 

A  s  I  — 


1  n  1 

- o  I  { $  -  2  T  a .  <|>  $  ( £ .  ) 

‘n= ![  n  (n  +  1 )]  2  ir.=-n  nm  i=l  1  nm  nm  1 


i  i 

+  y  y  a.a.<|>  ( £ .  )$  (  £  . )  i 


(3-26) 


Now  A  depends  on  I  parameters  {ai>  and  therefore,  the  op¬ 
timization  problem  requires  the  solution  of  a  linear  system  of 
dimension  I  +  1  : 


~  9(n  )*Q]p(n  )da(n  ) 


3 

3  X 


{ 


0 

(3-27) 
0  , 


which  we  will  write  in  the  form 


with  the  (1*1)  symmetric  positive  definite  Gram  matrix  Ut  , 

1  \47n^,  n  n  1  j  /  n  £  n  ( n +  1 )  J 2 

( 3-29a ) 


the  I-dimensi onal  vectors  U2>  x  ,  y] 


U2  =  [1,  1,. . .  ,  1}  . 


(3-29b) 


X  1  =  1  ’  a2 ’ ‘  ‘  ’  ai^  ’ 

yT  =  {  l  vi0)27rrr  2  *nm*m(Ci)} 

vn  =  1  m=-  n  ' 


and  the  scalars  x2  and  y2  , 


(3-29c) 


=  1 .  I  ,  ( 3-29d ) 


x2  =  A  ,  ( 3-29e ) 

y2  =  1  •  (3-29f) 

The  solution  of  the  system  (3-28)  can  be  easily  found  by  Cho- 
lesky's  method. 

3.2.  Approximation  of  degree  one  (J  =  1) 

The  next  simple  approximation  is  given  by  a  first  degree 
fit  which  requires  a  minimum  of  4  function  values  to  be  line¬ 
arly  combi  ned  .^However ,  let  us  consider  an  arbitrary  number 
2 

I  >  (C  +  1)  of  function  values;  then  Lf  will  be  approximated 
by 


1)  .. 


Lf  =  l  a  fU  )  . 

i  =  1 

.  =  the  dimension  of  the  space  including  zero  and  first 
degree  polynomials 


Then  Peano’s  theorem  states  that  the  remainder  Rf  can  be  re¬ 
presented  by 


Rf  =  jlC1(n)[<»*)(4*)|f(n)]do(n) 


(3-30 


with  (A*)  .  =  A*  +  2 
v  n 7  l  n 


(3-31 


and  the  spherical  Peano  kernel 


K,(n)  -  — R,G(l)U,n) 
1  (4«)2  S  ' 

and  Green's  function 


(3-32 


G j  ( 5  »n )  -  I  ni 1 

n  =  Z 


2n  +  l 


Pn(?-n)  • 


(3-33 


With  Schwarz's  inequality  we  obtain  the  estimate  for  1 Rf I 

'  'max. 

(Rf)2  s  J[K|(-i))2do(-<).|[(4*)(4*)|f(n))2do(n)  .  (3-34 

o  a 

As  in  section  3.1.  the  first  integral  will  be  denoted  by  A 
and  is  equal  to 


A  =  f[K  (n)]?do(n)  =  R,R  G(.2)U,c) 


/  1 .  (  4tt  )  4  < 

with  the  iterated  Green's  function 


(3-35 


Gf;U,n)  «  (4w)3  l 


2n  +  l 


n=2  [ n(n+l )]  [ n(n+l )-2] 


2  n 


P  U-n)  .  (3-36 


The  application  of  (3-35)  yields 


5  I  ([n(ntl)l[n(ntl)-21}'2  J  {  -  2*  J  a  *  U  ) 

n=2  m=-rr  i=l 

1  1  A 

+  y  y  a. a.*  (?.)♦  (5-)  \  (3-37) 

•  .  -  ,  i  jnminmj  I  v  ' 

i  =  I  .1  =  l 


The  factor  {[ n(n+l ) ][ n(n+l )-2] }”  goes  very  rapidly  to  zero 
and  as  a  consequence  only  very  few  terms  of  the  series  have  to 
be  considered.  Let  us  now  turn  to  the  evaluation  of  the  second 
integral  in  (3-34)  which  we  again  denote  by  B  , 


B  =  J((A*)(&*  +  2 )  f  (  n)  ]  2dcf{n) 


(3-38) 


With  the  representati on  (2-3)  for  our  kernel  we  obtain 
(a*)(a*  +  2)f  =  l  n ( n  +  1  )2[ n (n+1 )  -  2lanP(t) 

n=2 

and  considering  the  orthogonal i ty  relation  (3-15)  we  obtain 


b  -  4 ,i 


(3-39) 


The  relation  between  the  two  components  A  and  B  deserves  to 
be  discussed:  from  (3-10)'  and  (3-37)  we  conclude  that  the  number 
of  spherical  harmonics  and  its  local  integral  to  be  computed  de¬ 
creases  rapidly  with  increasing  degree  of  approximation;  at  the 
same  time  the  number  of  terms  needed  for  the  calculation  of  B 
increases  with  comparable  speed  (cp.  (3-16)  with  (3-39))  and 
with  (Rf)2  f  A*B  we  have  once  more  a  beautiful  example  of  the 
balance  of  difficulties:  blessing-burden  =  constant. 

The  best  estimation  of  the  coefficients  (a^  »  i=l»---»I 
is  also  getting  more  laborious  for  two  reasons:  fi rs t,  because 


of  the  relation  I  £  (0  +  1)  ,  and  second,  because  of  the 

2 

(J  +  1)  constraints,  the  linear  system  (3-28)  has,  for  the 

first  order  approximation  case  (J  =  1)  studied  in  this  section, 

2  2 

at  least  the  dimension  2(J  +1)  =  8  .  The  (J  +  1)  constraints 

are  due  to  the  required  annihilation  Rh  =  0  vh  e  ,  where 
P1  is  the  linear  space  of  all  zero  and  first  degree  polynomials 
in  three  variables,  restricted  to  a  .  Therefore,  the  optimiza¬ 
tion  problem  is  formulated  as  follows: 


TTF  ^  l  l  Xnm|jg(n)^nm(n)do(n)  -  I  a  ♦  (f  )~|}  =  0 

1  n=o  ra=-n  ^  D=  *  (3.40) 

lx  {  -  "  -  }  =  0 

nn  '  J 

which  leads  to  the  linear  system  (3-28)  with 

U ,=  L~l  v(,)P  ({..{.)}  ,  i,  j  -  1 . 1;  v(l):= - - J 

'  *?*n=2  n  "  1  3  J  "  [n(r.+l)l2[n(n+l)-2I2 

‘  ^  {4=  •  1  =  1-—  1  • 

^  (3-41) 


X|  =  {3^}  »  i  -  !*•••»  I  > 


y]  =  {  ^  vn)2r^T  l  *nini5i)}  ’  1=1 . 1  * 

vn=2  c=-n  * 

Note  that  x2  and  are  no  longer  scalars;  they  are  now  vectors 
of  length  4  ,  (because  of  the  4  constraints  in  the  linear 
model  we  need  4  Lagrange  multipliers). 
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yT  =  /4it(4> 


o  ,  o 


i  i  1  t  I  y 

1,-1  1,0 


As  in  section  3.1.  the  solution  can  be  obtained  by  applying 
Cholesky's  algorithm. 


3.3.  Approximation  of  arbitrary  degree  J 


The  general  case  of  approximation  degree  J  requires  a 

o 

linear  combination  of  I  a  (J  +  1)  function  values  and  the 
error  estimation  is  based  on  the  requi rement  that  the  function 
f  6  C2^J+,\o)  .  According  to  Peano's  theorem  the  remainder  Rf 
can  be  represented  by 


wi  th 


Rf 


o 


(3-42) 


(a  )  :  =.n  (a*  +  k  . )  =  (a*  +  <  )•••(&*  +  <  ), 
n 7 j  j=o'n  j  7  '  n  o7  v  n  j 


(3-43) 


and  the  eigenvalues  =  j(j  +  1)  of  the  Beltrami  differential 

operator  (Freeden,  1981,  p.  553).  The  spherical  Peano  kernel 
Kj(n)  is  given  by 


Ka(n) 


1 


(-4w) 

and  Green's  function  by 


j+T  V'»(c,n) 


(3-44) 


<\\  i  "  rJ  n  ~ 1 

Gj  (t.n)  =  (^)J  l  (Zn+1)  IT  (<-<•)  PnU*n)  (3-45) 

J  n=J+  1  U=o  n  n 


-  1 


(Freeden,  1981,  p.  556).  Evidently,  for  J  =  0  we  obtain  (3-7), 
for  J  =  1  equation  (3-33).  Schwarz's  inequality  provides  an 
estimate  for  |Rf|  , 

1  1  m  a  v 


(  Rf )  s  JtK^rOrdoU) 

o  i 

The  first  integral  is  equal  to 


[K)Jf'n)1  do(n) 


(3-46) 


A 


[Kj(n)] 


do(n) 


o 

with  the  iterated  Green's 


(4*r2(J+1)R  R  G®  u,o 

S  C 

function 


(3-47) 


GL2)(€.C)  =  (4tt)2^+1)  l  (2n+l)[~n  U  -<jl  2pn(^*4  ) » ( 3-48 ) 

n  =  J+ 1  n  n 

note  that  G*2)  is  simply  the  result  of  a  convolution  of  G(1) 

J 

with  itself. 


g(2)=  g(D»  g(D. 
j  j  j 


(3-49) 


With  (3-47)  and  (3-48)  A  can  be  shown  to  be  equal  to 
A  <  £  [".nfic  "ie.)l  y  (<j>2  -2$  y  a  *  (5  ) 

—  _  ,  b=0  n  L  \  nm  ynm,\  rnm'  1' 

n=J+ 1  ^  — 1  ra=-n  '  i=l 


+  l  l 

i=l  j=l 


a.a.<j,  U.)<J> 

x  3  nm  1  nm 


M  . 


(3-50) 


The  second  integral 
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reduces  to  a  very  simple  representation  if  f  is  a  homogeneous 
and  isotropic  function.  This  is  true  for  the  kernel  function 
discussed  here  (equation  (2-3)).  With 


(A*)  f(t)  -  f  1>~|(n*l)anPri(t)  (3-52) 


P  =  I  |[jno(Kn-<.)](n  +  l)ctn|2(2n+l)-1  .  (3-53) 


°»  r-  J 

l 

n  =  J  +  1 

and  the  orthogonal i ty  relation  (3-15)  we  obtain 

il.J 

n  =  J+  1 


The  best  estimates  of  the  coefficients  {a.}  are  obtained  from 

X 

the  solution  of  the  optimization  problem 


3a 


3  A 


ri 
{ 


J  n 


~  l  l  > 


n~o  ra--n 


nm 


g(n)-*.  (n)da(n)  -  l  a  ♦  U  ) 

J=1 


nm 

which  leads  to  the  linear  system  (3-28)  with 

i  “*  "2 

.n  (<  - < . 

h=n'  n  T 

n  =  J+  1 


}-° 

(3-54) 

}  =  o. 


Ui  {4*  l  GnoUn-^)  (2n  +  l)p  »  i  »j  =  !»•••»  1  » 

n=.1+1  l-  —  J  > 


U?  =  (c  )|  ,  n  -  0,...,  J, 

Z  ^  nm  1  J 


m  “n ) » « « n  )  i  !>•♦»}  I 


X  J  —  {  9  ^  }  y  i  If.'.  9  I 


C3-55) 


00 


l 

n  =  J+  1 


. n  ( k  -k  . 
J  =0  n  3 


l 

m  =  -ri 


1 


» 


» 


»  • 


y2  = 


n  =  0 , . .  .  ,  J  ,  m  =  -n,...n 


The  advantage  of  a  fast  converging  A  is  balanced  by  the  dis¬ 
advantages  of  a  slowly  converging  B  and  the  solution  of  a 
large  linear  system  Ux  =  y  .  Since  the  convergence  of  B  de¬ 
pends  strongly  on  the  value  of  a  and  therefore  on  the  depth 
D  ,  we  conclude  that  a  first  (probably  even  a  second)  degree 
approximation  is  advantageous  for  a  small  a  (large  D  ).  If 
a  is  very  close  to  1  (small  D  ) ,  we  are  obviously  in  troubles 
because  due  to  the  poor  convergence  of  B  we  have  to  keep  the 
degree  of  approximation  low  (probably  at  zero);  as  a  consequence, 
the  convergence  of  A  will  be  very  poor  and  the  estimation  of 
the  approximation  error  a  laborious  expensive  task. 


I 

t 

I 

* 

f 

!  4 .  GRAVITY  DISTURBANCE  MEAN  VALUES  VERSUS  DEPTH  OF  POINT 

MASSES 


In  this  chapter  we  shall  investigate  at  which  depth 
the  point  masses  are  most  likely  to  be  located  if  they  are  to 
be  derived  from  mean  gravity  disturbances  of  various  size  like 
5°  x  5°,  1°  x  1°,  etc.  In  other  words  we  look  for  (gravity 
disturbances)  averaging  operators  with  a  spectrum  (eigenvalues) 
as  close  as  possible  to  the  spectrum  of  the  operator  trans¬ 
forming  the  point  masses  into  gravity  disturbances. 

For  this  purpose  we  will  assume  a  continuous  (rather 
than  di screte )  anomalous  mass  distribution  on  a  geocentric 
sphere  with  radius  a  <  1  .  The  gravity  disturbances  are 

assumed  to  refer  to  the  geocentric  sphere  a  =  1  . 

According  to  equation  (2-3),  the  isotropic  integral 
kernel  K  of  the  operator  which  transforms  mass  anomalies 
into  gravity  disturbances,  can  be  represented  in  terms  of 

KU-n)  =  ^  l  (n  +  l)anP  U.n)  •  (4-1) 

ri=0  n 

We  employ  the  Funck-Hecke  theorem  (Muller,  1966),  which  states 

that  the  eigenvalues  k  of  an  isotropic  integral  kernel  K 

n 

on  the  unit  sphere  are  given  by 

l 

kn  =  2n  |  K( t)Pn(t)dt  (4-2) 

- 1 

with  P  denoting  the  Legendre  polynomial  of  degree  n  and 

n 

t  =  cosiji  .  Due  to  the  orthogonality  of  the  Legendre  polynomials 

l 

f  P  (t)P _<t)dt  =  JL 


6 


(4-3) 


N  S  *.-S 


the  eigenvalues  of  the  kernel  (4-1)  are  obtained  as 


(4-4) 


which  behaves  obviously  like  an  for  high  degree  n. 

The  question  is,  can  we  find  a  kind  of  moving  average 
operator  which  behaves  similarly? 

Let  us  investigate  the  moving  average  operator  acting 


on  a  circular  cap  of  radius  4,  , 

HO  =  {  BU.n)f(n)do(n-) 


with  the  isotropic  integral  kernel 


(4-5) 


B(c • n  )  = 


1  for  £.n  ^  cos 


0  else 


(4-6) 


According  to  the  above  referenced  Funck-Hecke  formula,  the 
eigenvalues  &n  of  B  are  given  by 


B„  ■  T^T  f  p„^)dt  •  <4-7> 

o 

to 

The  integral  yields  Tfffj  lPnM(tQ)  ”  ^n+1^0^  *  which  can 

be  used  to  design  a  recursion  formula  for  8  (Sjoberg,  1980) 

n 


8l(to>  =  7  (1  +  to> 


(4-8) 


»„<*<>> 


(n-2 ) B  (t  )]  ,  n  s  2 

n- 1  o 


These  eigenvalues  approach  zero  by  oscillating  around  zero, 
dependent  on  tQ  .  For  the  extreme  case  ^ o  =  0  (no  smoothing), 
all  eigenvalues  are  equal  to  1  ,  for  the  extreme  4>q  =  * 

(total  smoothing),  all  eigenvalues  are  zero  apart  from  the 
zero  degree  eigenvalue  which  is  equal  to  1  . 

Now,  the  behavior  of  these  eigenvalues  does  not  quite 
agree  with  that  of  kn  given  by  (4-4);  an  approximation  pro¬ 
cedure  could  be  designed  which  fits  (4-4)  to  (4-8)  in  the  best 
possible  way  yielding  a  best  possible  value  for  a  and,  there¬ 
fore,  for  the  depth  of  the  mass  anomaly  layer.  For  the  purpose 
of  obtaining  a  closed  expression  for  the  covariance  function 
of  mean  values  such  an  approximation  is  actually  used;  this 
is  basically  equivalent  to  replacing  the  mean  values  at  zero 


level  by  point  values  at  a  certain  altitude  dependent  on  a  . 
This  technique  provides  us  with  information  on  how  to  select 
the  depth  of  a  mass  layer  corresponding  to  gravity  disturbance 
mean  values  of  rectangular  blocks  of  a  certain  size: 


a 


block  size 


depth 


Denoting  the  harmonic  coefficients  of  the  mass-distur¬ 
bance  implied  gravity  disturbance  6g  by  g  and  the  har- 

nm 

monic  coefficients  of  the  distribution  of  mass  disturbance 

2irip  by  p  ,  we  conclude  from  (4-4)  that  they  are  related 
nm 

to  each  other  by 


nm 


n+1 

n 


n 

a  u 


nm 


(4-9) 


Since  a  <  1  ,  the  behavior  of  the  mass-disturbance  implied 
gravity  disturbance  at  zero  level  will  be  smoother  than  that 
of  the  mass  disturbances  with  the  degree  of  smoothness  depen- 
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ding  on  the  depth  of  the  mass-layer.  Vice  versa,  deriving 
mass  disturbances  from  unsmoothed  gravity  disturbances  is 
an  unsmoothing  process  and  comparable  to  differentiation. 

A  stable  transformation  from  one  quantity  to  another  one  can 
be  expected  if  the  frequency  behaviour  of  the  output  function 
is  smoother  than  or  at  least  as  smooth  as  that  of  the  input 
function.  Therefore  it  is  a  good  advice  to  use  mean  gravity 
disturbances  6g  with  harmonic  coefficients  *5  related  to 

no 

9nm  by 


6 

n 


nm 


(4-10) 


to  determine  6 u  such  that  the  behaviour  of  k^  is  comparable  to 
that  of  (4-8).  In  other  words,  the  unsmoothing  6g  •+  6p  has  to  be 
balanced  by  a  sufficiently  smooth  6$. 


As  far  as  the  determination  of  a  "best"  value  for  the  factor 
ij>o  is  concerned,  one  could  follow  the  ideas  of  Schwarz  (  1976) 
and  fit  the  eigenvalues  {  kn>  best  (in  the  sense  of  1 eas t  -  squares 
e.g.)  to  the  eigenvalues  {  Bn}  with  n  _<  N  (a).  If  a  belongs  to 
a  large  D,  N(a)  will  be  low  (on  the  order  of  a  few  hundred  for 
0^500  km,  if  a  belongs  to  a  small  D,  N(a)  will  be  high  (on  the 
order  of  a  few  thousand  for  D  'vlO  km);  the  degree  of  truncation 
N ( a )  depends  on  the  depth  of  the  point  mass  layer  D.  (See 
chapter  2  of  SUnkel  (1981).) 
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In  (Sunkel ,  1982)  it  was  argued  that  a  mass  distribution 
confined  to  a  single  layer  is  inadequate  for  the  representation 
of  the  disturbing  potential.  Mass  distributions  on  a  few  layers 
at  various  depths  a  re  pre f erabl e .  In  chapter  4  we  have  discussed 
how  to  determine  the  most  appropriate  depths  provided  a  cer¬ 
tain  mean  value  information  on  gravity  disturbances  at  zero 
level  is  available. 

Let  us  now  describe  the  principle  of  the  procedure  which 
could  be  used  to  derive  point  masses  from  known  gravity  dis¬ 
turbances.  For  the  sake  of  simplicity  we  will  assume  that  the 
various  mean  values  can  be  approximated  sufficiently  well  by 
a  corresponding  moving  average  as  described  in  the  previous 
chapter.  Then  the  actual  gravity  disturbance  harmonic  coeffi¬ 
cients  (unsmoothed)  can  be  represented  in  terms  of 


=  [l  <e< 

Li=o 


(1)  aU+l) 
n  "6n 


(5-1) 


with  e(o)  =  1  vneH  and  b(l+1)=  6  .  Let  us  illustrate 

n  o  n  n  ,  o 

equation  (5-1)  by  an  example:  consider  gl  '  as  the  eigenvalues 

n 

of  the  moving  average  operator  correspondi ng  to  a  mean  value 


of  10  '  x  10'  ,  b 


correspondi ng  to  1°  x  1°  ,  and  g 


correspondi ng  to  5°  x  5°  .  Then  ( g (L) -6 ( L+ 1 ) )g  =  g{3)g 

n  n  nm  n  nn 

would  be  the  coefficients  of  the  5°  x  5°  mean  values, 
(g(2)-g(3))9  the  coefficients  of  the  1°  x  1°  mean  values 

n  n  nm 

referred  to  the  5°  x  5°  mean  values,  )g  the 

n  n  nm 

coefficients  of  the  10'  x  10'  mean  values  referred  to  the 
1°  x  1°  mean  values,  and  (glo)-e(1>)9  =  (1-  sj^g  the  co- 

efficients  of  the  actual  gravity  disturbance  field  referred 
to  the  10‘  x  10'  mean  values. 


V 


r 


According  to  (.4-9)  the  coefficients  g 

nm 

coefficients  of  the  mass  disturbances  by 


are  related  to  the 


9 


nm 


n  +  1 


(1 ) 
rim 


(5-2) 


The  coefficients  y(1)  refer  to  the  mass  disturbances  at  the 

n  m 

layer  1  ;  the  layer  1=0  is  assumed  to  be  the  most  shallow 
and  1  =  L  the  deepest  layer.  Comparing  (5-2)  to  (5-1), 


l  (6. 


(1+1) 


l  n  +  ‘ 

1=0  n  +■> 


n  (1) 
aiMnm 


(5-3) 


then  a  "separation  by  layer  and  mean  value"  presents  itself 
as  a  practicable  method  to  determine  mass  disturbances  from 
mean  gravity  disturbances. 


(0<1)_0C1+1> 

n  n 


= 


n  Cl) 
a ,  y 
1  nm 


(5-4) 


In  other  words,  we  intend  to  determine  the  mass  disturbances 
located  at  the  layer  1  from  gravity  disturbance  mean  values 
corresponding  to  the  smoothing  1  and  referred  to  the  (1+1) 
smoothed  values.  For  our  previous  example  this  means  that  the 
point  masses  at  depth  0  =  D3  (1  =  L  =  3)  should  be  de¬ 

termined  from  5°  x  5°  means,  the  point  masses  at  depth 
D  »  D£  (1=2)  from  residual  1°  x  1°  means  (referred 

to  5°  x  5°  means),  the  point  masses  at  depth  D  =  Dj 
(1  =  1)  from  residual  10'  x  10'  means  (referred  to  1°  x  1° 
means),  and  the  point  masses  at  zero  depth  D  =0  from  the 
residual  gravity  disturbances  (referred  to  the  10'  x  10’ 
mean  values). 
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s-1-  Po-1nt  masses  of  50  X  5°  and  1°  X  1°.  etc .  global  mp*n 
gravity  dis turbances 


It  is  natural  to  put  mass  points  at  the  center  of  each  mean  value 
block.  Therefore,  to  each  5°  x  5°  mean  value  there  corresponds  one 
mass  point  at  depth  0  =  DL;  given  the  vector  of  5°  x  5°  mean  gravity 
disturbances  g,  we  have  to  find  the  vector  of  point  masses  p  by 
solving  the  linear  system  of  equations 


g  =  Cp 


(5-5) 


Let  us  now  investigate  the  structure  of  the  matrix  C 

in  detail.  We  assume  that  the  sphere  is  subdivided  into  I 

(equally  spaced)  parallels  and  J  equally  spaced  meridians 

and  that  to  each  grid  element  [$. ,  ;  x.  ,  A.  ]  there 

1  1+1  3  j  +1 

corresponds  one  mean  gravity  disturbance  to  its  center  of 
point  mass.  Then  g  and  p  consist  of  I  subvectors  (cor¬ 
responding  to  I  parallels)  consisting  of  J  mean  values  and 
point  masses,  respectively,  (corresponding  to  J  meridians); 
in  the  same  way  the  matrix  C  consists  of  I2  submatrices 
of  dimension  J2  each,  correspondi ng  to  the  subdivision  of 
the  data  vector  g  and  the  solution  vector  p  , 


g  =  fg 


t11  C12 

c2i  c22 


C11  (5-6) 


.ii  ' 


C11  is  the  matrix  which  relates  g1  to 
coefficient  which  relates  gM  to 
distribution  of  the  data  and  the  unknowns  with  respect  to 


v~  ,  c*r,  is  the 

j  j 

From  the  regular 


longitude  it  is  evident  that  c 
tude  difference  a 


1 1 
3  j  ‘ 


depends  only  on  the  longi- 
-  A.  for  each  pair  of  parallels  ($. 
Matrices  of  this  kind  are  known  as  Toeplitz  circulant  or 
briefly  circular  matrices,  because  row  number  j  equals  row 
number  j-1  rotated  by  1  element.  Circular  matrices  become 
diagonal  matrices  under  a  discrete  Fourier  transformation. 

This  property  can  be  used  very  advantageously  to  design  a  very 
efficient  algorithm  for  the  solution  of  the  very  large  linear 
system  (5-5),  using  frequency  domain  methods.  In  the  sequel, 
we  shall  briefly  outline  this  method  following  closely  the 
fundamental  work  by  Colombo  (1979). 

The  first  row  (and  all  the  others)  of  every  submatrix 
is  an  equispaced  sample  of  an  even  function;  therefore,  it  has 
only  a  non-vanishing  cosine-spectrum  with  Fourier-coef f icients 
c^1*  (an  overbar  denotes  a  quantity  in  the  spectral  domain, 
k  stands  for  a  discrete  frequency). 


*C  i  i  '  * 


ff 


J-1 

l 

\  ’  =0 


.  1 1 
'°j 


c“, cosuik j  ’ 


(5-7a) 


and  vice  versa,  the  elements  of  this  particular  row  are  obtained 
by  an  inverse  Fourier  transform. 


r^- 


-  30 


c11 ,  =  l  c^1  coswkj 1 
03  k=o  k 

with  w:  -  2it/J  , 

K:  =  i n  t( J/2 )  , 


C5-7b) 


( 5-7c ) 


=  -l 


if  k  =  0  or  k  =  K  (J  even) 


^  J/2  if  0  <  k  <  K  ( =K  if  J  odd)  . 

The  row  number  j  equals  row  number  j-1  ,  rotated  by  one 

element,  and  therefore  ci£ ’  i s  synthesi zed  by 

j  j  ’ 


=  l  c*1 '  coso>k(  j  '  -j  ) 
33  k.=o 


(5-8) 


Using  the  orthogona 1 i ty  relations  between  equispaced  sampled 
trigonometric  functions 

0  if  k  =  k'  =  0  or  k  =  k'=K  (J  even) 

£  cosukj  cosuk'j  =  ■  J/2  if  0  <  k  =  k'  <  K  (=K  if  J  odd) 

3-°  0  if  k  /  k' 


j_i  J/2  if  0  <  k  =  k'  <  K  (=K  if  J  odd) 

l  sinwkj  sinu  k'j  =  ■ 
j=o  0  else 


1  cosukj  si  no,  kj  '  =0  v  k,  k'  , 

j=o 


(5-9) 


we  obtain 


, , .  cosukj  j-i  cosukj ' 

c“  H  =  l  Cxi  \  | 

k  sinukj  j ’ =o  33  sinwkj’ 

l  J  l  J 


,  0  s  k  ^  K  .(5-10) 


Consequently 


c*  :  =  {coswkj} 


s£  :  =  {sinwkj} 


i  j  ■  0,  .  •  .  ,  J— 1 


(5-11 


are  eigenvectors  and  c^  H  eigenvalues  of  the  submatrix  C; 

Let  us  now  transform  the  vectors  g  and  p  into  the 
frequency  domain  ( gf,  7)  » 


9k  =  <  • 


k  k 


with  g1  a  yi|  k 
k  k  s 

k 


-i  _  ,i  k 

V  '  V  s  , 

k 


(5-12 


(5-13 


observing  (5-5)  and  (5-10),  we  find  immediately  the  relation 

between  v1  and  61  , 
k  k 


1  =  y 

k  i  •  =  1  k  k 


(5-14 


Denoting  the  1*1  matrix  of  eigenvalues  cL1  H  by  E  , 

k  k 

equation  (5-14)  is  given  by 


Here  we  have  already  the  algorithm  at  hand  which  permits  the 
fast  transformation  from  the  known  mean  gravity  disturbance 
vector  g  to  the  unknown  point  masses  v  : 

1)  Calculate  the  first  row  of  each  submatrix  CiL' 
and  transform  it  into  the  frequency  domain  by  fast 
Fourier  transformation; 

2)  Establish  the  matrices  E  for  k  =  0,  ....  K 

k 

and  calculate  the  inverses; 

3)  transform  the  data  vector  g  into  the  frequency 
domain  (y  )  i 

*  1 

4)  perform  the  matri x-vector  product  6  =  E  y 

k  k  k 

for  k  =  0 ,  .  .  .  ,  K  ; 

5)  Perform  the  inverse  Fourier  transformation  of  & 

k 

(synthesis)  and  obtain  p  ,  the  vector  of  point  masses, 

M1  -  l  K  -  (5—16 ) 

k=o 

The  advantages  of  using  this  frequency  domain  method 
are  obvi ous : 

(a)  Instead  of  solving  the  linear  system  (5-5)  with  dimen¬ 
sion  I J - 1 J  ,  we  have  to  solve  the  K+l  linear  systems  (5-15) 
of  dimension  I»I  each.  Since  the  solution  time  is  approxi¬ 
mately  proportional  to  the  3rd  power  of  the  dimension  of  the 
linear  system,  the  computation  time  reduces  by  a  factor  of 
about  J2  ( -5200  for  the  5°  x  5°  case,  -130  000  for  the 

1°  x  1°  case). 

(b)  Since  all  subma trices  C11’  are  Toeplitz  circulant, 
we  need  only  calculate  the  first  row,  and  moreover,  since  the 
first  row  represents  an  equispaced  sample  of  an  even  function, 
only  half  of  the  elements  have  to  be  calculated.  This  reduces 
the  (anyway  expensive)  calculation  effort  for  the  setup  of  C 


33 


by  another  factor  2J  compared  to  the  straight-forward  al- 
gori thm. 

(c)  The  storage  requirements  are  drastically  reduced  due 
to  redundancy  in  the  C-matrix  . 

(d)  The  calculation  of  all  required  elements  requires  the 
evaluation  of  an  integral.  In  chapter  3  we  have  proposed 
optimal  approximation  algorithms  for  its  calculation,  which 
can  be  easily  implemented. 

Note  that  in  the  spherical  case  C  is  not  symmetric 

£  £  *  £4  • 

due  to  the  convergence  of  meridians  (c....,  t  c....  because 
the  coefficient  which  relates  a  mean  value  in  row  i  to  a 
point  mass  in  row  i'  /  i  is  different  from  the  coefficient 
relating  a  point  mass  in  row  i  to  a  mean  value  in  row  i’  ; 
for  a  cylinder,  C  would  be  symmetric.)  However,  if  the 
arrangement  of  mean  values  and  point  masses  is  symmetric  with 
respect  to  the  equator,  the  number  of  different  elements  of 
C  is  once  more  reduced  by  a  factor  of  2  . 
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